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Abstract 

We propose time-domain boundary integral and coupled boundary integral and 
variational formulations for acoustic scattering by linearly elastic obstacles. Well 
posedness along with stability and error bounds with explicit time dependence are 
established. Full discretization is achieved coupling boundary and finite elements; 
Convolution Quadrature is used for time evolution in the pure BIE formulation and 
combined with time stepping in the coupled BEM/FEM scenario. Second order 
convergence in time is proven for BDE2-CQ and numerical experiments are provided 
for both BDE2 and Trapezoidal Rule CQ showing second order behavior for the 
latter as well. 
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1 Introduction 

The study of the acoustic scattering by a linearly elastic obstacle and its corresponding 
elastic response has been a subject of interest in both the mathematical and engineering 
community for some years now. In the case of time-harmonic regime, the study of the 
existence and uniqueness of solutions dates back at least to 1986 [21]. The well-posedness 
of several purely boundary integral formulations was analized in [251123 El] by assuming 
that the scatterer had at least a boundary of class C^. Later on, in the 2000’s, combined 
boundary integral and variational formulations were proposed in [26] |T2] and proven to 
be well posed also for smooth scatterers. In these works the elastic response is modeled 
variationally and the unbounded acoustic scattering is treated with either a boundary 
integral equation or by introducing an artificial boundary where an absorbing condition 
is imposed. Recently, the more general case of a Lipschitz scatterer was dealt with in [7] 
within the framework of a variational formulation with a hctitious boundary. 
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On the computational side, the coupling of Boundary Elements and hp-Finite Elements 
was studied in [8] where the Burton-Miller equation is used to model the acoustic waveheld; 
the authors provide a posteriori error bounds aiming for an adaptive implementation. 
The hctitious boundary approach with hnite elements has been thoroughly investigated 
in [33l Uni EO] and a DG-like implementation was carried out recently in |6]. 

The transient regime, on the other hand, has not enjoyed so much attention -at 
least in the mathematical community- as its frequency-domain cousin. In [161 E] fhe 
problem is posed in a slab-like unbounded domain imposing Erst order absorbing bound¬ 
ary conditions, while in [2Z] well-posedness is established for the coupled boundary inte¬ 
gral/variational formulation also in a slab-like region. Within the engineering community, 
the time-domain case has attracted attention at least since 1991. Numerous approaches 
have been attempted without much theoretical justihcation but with satisfactory results. 
To cite some examples, BE/FE coupling with Convolution Quadrature was employed 
in [13], BE/BE coupling using Newmark time integration was the preferred treatment 
in [22] and FE/FE coupling with an absorbing boundary condition and Newmark time 
integration were used in [TH|. A comprehensive list of related work can be found in [37] . 

The present work strives to £11 the gap in the mathematical analysis of the time-domain 
wave-structure problem. It deals both with the pure boundary integral formulation -which 
arises naturally when dealing with homogeneous acoustic and elastic domains- and the 
combined boundary integral/variational treatment where integral equations are used only 
for the acoustic dynamics, being better suited for general elastic scatterers. The former 
case leads to a numerical treatment where only Boundary Elements are used for space 
discretizations while the latter is naturally suited for a coupled Boundary Element/Finite 
Element implementation. 

Despite the fact that each formulation requires a very different numerical discretiza¬ 
tion, the techniques and tools required to carry out the theoretical study are surprisingly 
similar. Following [23 ESI E] the analysis is done in the Laplace-domain aiming for a 
Convolution Quadrature treatment similar to that done for the purely acoustic case in 
[3 [HE]. We deal simultaneously with the continuous and discrete cases by posing the 
problems in a general closed subspace of the appropriate function spaces. Well-posedness 
is proved in the Laplace-domain via an equivalent exotic transmission problem for which 
a variational formulation is found. The resulting stability bounds are written carefully in 
terms of the Laplace parameter s in order to apply results from [35] which give explicit 
time-domain estimates. Error bounds in the time-domain are obtained following a similar 
approach for the semidiscrete problem. Full discretization and convergence estimates are 
given for the case of BDF2-CQ. 

Finally, numerical experiments are carried out using Boundary Elements coupled with 
BDF2 and Trapezoidal Rule Convolution Quadrature for the BIE formulation and BE/FE 
with Trapezoidal Rule Convolution Quadrature for the boundary element part paired 
with Trapezoidal Rule time stepping for the elastic field. The results support the order 
of convergence predicted in the theoretical part of the paper. 

A word on real and complex Sobolev spaces. Basic knowledge of Sobolev spaces 
its trace space and its dual is assumed throughout. In 

everything that follows, the Sobolev spaces H^{Vt) and will be used with the 
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same notation for real-valued and complex-valued functions. If, for instance, we take 
a subspace V C we will understand that it is a subspace of the real-valued 

space and that its complexification will be used whenever complex values are 

considered. Scripted parentheses will be used for real inner products of scalar, vector- 
or matrix-valued functions: 

(a,6)n:= ab, (a, b)n:= / a ■ b, (A,B)q:= / A : B, 

Jfi Jfi Jq 

where in the latter the colon denotes the Frobenius inner product of matrices. When the 
functions take complex values, we will still use the parenthesis with this precise meaning 
and will explicitly conjugate quantities whenever needed. When we change the font de¬ 
scribing a Sobolev or Lebesgue space from italic to boldface we mean the product space 
of d copies of the space. For example, L^(r2) := 

2 Homogeneous isotropic solids: BIE formulation 

The PDE system. Let hi- C be a bounded, not necessarily connected region, 
lying on one side of its Lipschitz boundary F, and let 17+ := \ hi. be its unbounded 

complement. The unit normal vector held on F, exterior to 17-, will be denoted u. Our 
problem can be explained as follows: an incident acoustic held arrives at an obstacle 
at time 7 = 0 and interacts with a homogeneous isotropic elastic body occupying f7_. The 
elastic properties of the material are represented by the two Lame parameters that dehne 
the structure of the linear stress tensor: 

cr(u) := 2/i£(u) -h A(V ■ u)I, £(u) := |(Vu -F (Vu)'*'), 

where I is the d x d identity matrix and u : f7_ —)• is the displacement held. This will 
be the only occurrence of the Lame parameters (until we get to the section on numerical 
experiments), and we will use the Greek letters A, /i to represent other quantities. Related 
to the stress tensor we can dehne the Lame operator A*u := V ■ cr(u) and the normal 
traction on the boundary by t(u) := cr{u)ix. Before we give a rigorous mathematical 
formulation of the problem, let us start by writing the system of PDE with transmission 
conditions that we want to solve: 


psWt = A*u 

in f7_ X [0, oo). 

(2.1a) 

c~‘^vtt = An 

in 17+ X [0, oo). 

(2.1b) 

■Ufzz=a,(n + n‘“) 

on F X [0, oo). 

(2.1c) 

t(u) = -p/(n-t-n“'=)fjz 

on F X [0, oo). 

(2. Id) 


Here pf and ps are the respective constant densities of the huid and elastic media, the t 
subscript denotes partial diherentiation with respect to time, and dy is the normal deriva¬ 
tive operator on F. This system is complemented with homogeneous initial conditions for 
u and V (and their time derivatives), and a causality condition that can be expressed as: 
for all 7 > 0, n = 0 except in a bounded region (that changes with 7). A derivation of this 
model can be found in [2H]- 
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Figure 1: A cartoon of the geometric setting: the solid (brown) is surrounded by the 
unbounded medium (blue) 


Traces on the boundary. In order to properly dehne our concept of solution to (2.1), 
we will need some additional notation. Given u G H^(M'^ \ T), we consider its interior, 
exterior, averaged, and difference traces: 

7"u, 7 +u, ^7u}} := |(7“u + 7+u), M := 7 “u - 7 +u. 

For u G H^(]R‘^ \ F) such that cr(u) G L‘^{W^ \ F)'^^'^ we dehne the weak interior and 
exterior traction helds using Betti’s formula 


(FF(u),7v)r := ±(cr(u),£( v))q^ ± (A*u, v)q^ Vv G H^(M‘^). 

These are elements of the dual space H~^/^(F). From now on the F-subscripted an¬ 
gled bracket will be used to denote the duality product of with Ff^/^(F) or 

H-V2(r) 

and H^/^(F). We will also dehne the average and jump of the normal traction 
as ■§)t(u)(§- := |(t+(u) -|-t“(u)) and |t(u)] := t~(u) — t+(u). Two sided traces (and 
normal derivatives) for scalar functions (with Laplacian in L^) are similarly dehned. The 
following two operators related to the normal vector held 

N : Hi/2^F) ^ N* : ^ H^/^^F) 

(p I —> (p ■ i', (p I —> (pi', 


will be used to give rigorous meaning to some elements appearing in the transmission 
conditions. 


Weak form. For the sake of completeness we will now give the weak form of the equa¬ 
tions (2.1). Note that all the estimates that we will produce will be developed using the 
Laplace transformed equations, and it will be only those equations that we will need to 
deal with rigorously. We look for a pair (u, v) of causal distributions, with values in the 
space 

{u G Hi(G_) : A*u G L2(G_)} x {n G H\n+) : Av G L‘^{n+)} 

such that 


Psii = A*u 

(in L2(G_)), 

(2.2a) 

c~‘^v = An 

(in L\n^)), 

(2.2b) 

7 “u ■ u = dp^v + tto 

(in 

(2.2c) 

f(u) = - pfi'y+i) + l3o)i' 

(in H^/\T)). 

(2.2d) 
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In (2.2), the upper dots are used for distributional differentiation and the parentheses in 
the right-hand sides tell where the distributions are compared. Full details on how to 
understand wave equations in the sense of vector-valued distributions can be found in 
[53] . Also, we have used cto and /3o to denote general causal distributions with values 
in and respectively. Existence and uniqueness of solution to (2.2) can 

be proved with some additional constraints: we have to assume that the data and the 
solution are Laplace transformable with Laplace transforms dehned in a subset of the 
form Res > (Jo for some (Jq. 


Laplace-transformed system. Let us now consider a slightly different problem. Now 
Ao G and (j)o G are data, and we look for (u, n) G x 

such that 


2 A * 

Py:S u = a u 
{s/cf‘v = An 
S7“u ■ u = d^v + Ao 

t~(u) = - p/s(7+n + (l)o)h> 


in f2_. 

(2.3a) 

in f2_|_. 

(2.3b) 

in F, 

(2.3c) 

in F. 

(2.3d) 


This problem will be studied for all s G C+ := {s G C : Res > 0}. The relation between 
(2.3) and (2.2) is simple: if Ao = £{q;o}(s) and 0o = 'C{/5o}(s), then the solution of (2.3) 
is the Laplace transform of the solution of (2.2). 


Calderon calculus for the acoustic problem. The single and double layer operators 
associated to the Laplace resolvent equation (the Laplace transform of the wave equation) 
can be dehned as the solution of a transmission problem. For given ( 75 , 7 ) G x 

^-i/2(f) g g problem of hnding v G \ F) satisfying 

An — (s/c)^n = 0 in \ F, 

[ynl = 7?, 

{d^vj = r], 

has a unique solution, which we write in terms of two linear operators, known as the single 
(S) and double (D) layer potentials 

n = S(s/c )7 — D(s/c)75. 

Associated to the potentials, we can dehne four integral operators 

V(s) := ^7 ■ }}S(s) = 7S(s), K(s) := §7 ' SD(s), 

K*(s) := id, ■ SS(s), W(s) := -id, • }}D(s) = -d,D(s). 


Calderou calculus for the elastic problem. In parallel to the dehnitions for the 
acoustic case, we can dehne the layer potentials and operators for the resolvent Navier- 
Lame equations (Laplace transforms of the linear elasticity equations) by solving a trans¬ 
mission problem. For given G H^/^(F) x H~^/^(F) and any s G C+, we hnd 
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u G \ r) satisfying 


A*u — = 0 in \ r, 

[yu] = ip, 

[t(u)l = T], 

and write the solution in terms of the data 

u = S(s)r/ — D(s)(/? 

by introducing the elastic layer potentials. Finally, we can dehne four integral operators 

V(s) := {{7 ■ SS(s) = 7 S(s), K{s) := §7 ' Ws), 

K*(s) := iti ■ )}}S(s), W(5) := -iti ■ )SD(.) = -tD(5). 


The boundary integral system. The boundary integral system equivalent to (2.3) 
is derived by choosing 0s := 7 ~u and <pf '■= as unknowns, using the representation 
formulas for v and u and hnally imposing the transmission conditions. The process is 
quite standard and we will only sketch the main steps. We introduce the matrices of 
operators 


L(s) 


W(s) + p/s2NW(s/c)N p/s(N'K(s/c) - K*(s)N') 
Pfs(NK{s) - K*(s/c)N) (pjs)2NV(s)N' + PfW{s/c) 


and 


M(s) 


-p/sNW(s/c) p/s(-|l + K*(s))N* 
Pf{ll + K\s/c)) -(p/s)2NV(s)N* 


Denoting ]HI=‘=^/^(r) := x it is easy to show that by well-known prop¬ 

erties of the boundary integral operators on Lipschitz domains (see the general theory in 
[31]) L(s) : eV2(r) ^ e-i/2(r) and M(s) : iWV2(r) X ^ e-h2(r) are bounded. 

For the sake of notational simplicity, we will write L(s)(0, 0), understanding that the 
vector (0, 0) is hrst transformed into a column vector and them left-multiplied by L(s). 


Theorem 2.1. //(u,n) solves (2.3), then (0e,0/) := (7 u, 7 +n) satisfies 


L(s)(0s,0/) = ]R(s)(Ao,0o). 


Reciprocally, if fif) is a solution of (2.4), then 


u := - p/sS(s)N*(0o -t- 4>f) - D(s)0s, 
V := S(s/c)(Ao -h sN0s) -h D(s/c)0/, 


(2.4) 


(2.5a) 

(2.5b) 


is a solution of (B- 

Proof. If (u, n) satishes ( |2.3a[ )-( p.3b[ ), then we have the representation formulas 

u = S(s)t~(u) — D(s) 7 “u and v = —S{s/c)dfv + D{s/c)''y~^v, (2.6) 
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and the boundary integral identities 


^7 u = V(s)t (u) — K(s )7 u, |t (u) = K*(s)t (u)+W(s )7 u, (2.7a) 

17 + 1 ; = —V(s)c?j'n + K(s) 7 ’''n, \d^v = —K*(s)c?J'n — W(s) 7 ^n. (2.7b) 


If we define (0s, 0/) := (7 u, 7 +n) the transmission conditions (2.3c)- (2.3d) become 

t~(u) =-psN*(0/-1-0o) and =-(sN0s + Aq). ( 2 . 8 ) 


Substituting (2.8) in (2.7) and adding the two equations in (2.7a) and the two in (2.7b) 


gives the integral equations (2.4) 


Reciprocally, let (u,n) be dehned by (2.5) where (0s,0/) solve (2.4). Since (u,n) are 
dehned with potentials, it follows that (2.3a) and ( |2.3b ) are satished. Applying the limit 
formulas (2.7|) in (2.5), we see that 


f (u) p/sN‘( 7 +n + 0o) = - fpjsN‘(0o + 0/) - P/sK*(s)N*(0o -7 0/) + W(s)0s 

-h p/sN*V(s/c)(Ao -h sN0s) p/sN*(|0/ K(s/c)0/) 

= (W(s) + p/s2nV(s/c)N)0s + p/s(N*K(s/c) - K*(s)N*)0/ 
+ pfsWY{s/c)\o + P/s(iNVo - K*(s)NVo) 

= 0 , 


by the hrst equation in (2.4). This proves the hrst transmission condition (2.3c). Similarly 
(2.3d) is proved using the second equation in (2.4). □ 


3 Numerical discretization of the BIE system 

3.1 Stability of Galerkin semidiscretizations in space 

Galerkin semidiscretization in space. We next consider a Galerkin discretization 
of the integral equations ( |2.4 ). Note that when returning to the time-domain (by taking 
inverse Laplace transforms) this is simply a Galerkin semidiscretization in space of the 
system of delayed integral equations whose Laplace transform is (2.4). Following |29], the 


study of solvability for (2.4) is done at the same time as the study of Galerkin stability. 


We then choose two closed subspaces Yh C Hi/2(r) and W C For Galerkin 

semidiscretization, these spaces are taken to be hnite dimensional. In the case of non¬ 
discretization (analysis of well-posedness of (|2.4|)) we just take Y^, = and 10 = 


The Galerkin discretization of (2.4) seeks (0|,0j) G Y^ x 10 satisfying 


(L(s)(</>|,7),(A<\7)>r = (R(s)(A„,0„),(^\7))r V(m\7) € x Y„. (3.1) 

The angled bracket is the duality product of ]HI“^/^(F) with Wg gan also write 


(3.1) in the very condensed form 


L(s)(,/.|, ,1,'}) - E(s)(A„, ,/>„) e Y; X = (Yi X Yi,r, 


(3,2) 
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where X° denotes the polar set of X, that is, the set of elements of the dual space that 
vanish on X. Following the same techniques of [29] we hrst rewrite the Galerkin equations 


(3.1) as an exotic transmission problem. Note that in the new transmission problem, the 
elastic and the acoustic fields are dehned on both sides of the boundary. 

Proposition 3.1 (Transmission problem for Galerkin equations). Let {4>^, cf/j) EYhXYh 
satisfy (3.1) and let 



u* := -p^sSWN‘(.#,(, + 7) - D(s)</.|, 

(3.3a) 


:= S(s/c)(Ao + sN0|) + D(s/c)0j. 

(3.3b) 

Then (u^, 

G \ F) X \ F) satisfies: 



= 0 mM‘^\F, 

(3.4a) 


Xv^ — [s/cYv^ = 0 m \ F, 

(3.4b) 


sN[7ui - = -Ao, 

(3.4c) 


[t(u^)l - p/sN^Iyn'^] = -p/sNVo, 

(3.4d) 


(hu'-|,|7i,'>i) e Y^xn, 

(3.4e) 


(sNy+u^ + a-n^,t+(u'*) + pfsW-f-v^) G Y° x Yf. 

(3.4t) 

Reciprocally, 

given a solution of (3.4) the functions 



(</.|,y) ^^Ihul.-ht’DeYixyi. 

(3.5) 


satisfy (3.1). 


Proof. It is clear that the functions dehned by (3.3) satisfy (3.4a) and (3.4b). Moreover, 
lyu^] = 02 and [[yn^] = —0/, and therefore (3.4e) is satished. At the same time 

|t(u'^)l = -p/sN‘(0o + 0j) = -p/sN*(0o - 
{dyV^l = Ao + sN 02 = Ao + sNlyu'^], 

which proves ( 3.4c| ) and (3.4d). Finally, using the jump properties of the potentials, it is 
easy to verify that 

(sN 7 +u'^ + a^n^,t+(u'*)+p/sN* 7 "n'^) = L(s)(02, 0/) - I^(s)(Ao, 0o), (3.6) 

which proves ( |3.4f| ) (see ([T^). 


Reciprocally, if we are given a solution of (3.4) and we dehne (0|, 0j) with (3.5), then 
by the representation formulas and (3.4c)-(3.4d), it follows that we can write the helds 
(u^,n^) as in (3.3). We can then use (3.6) again and prove that (3.4f) implies (3.1). □ 


The next step consists of hnding a variational formulation for (3.4). This will be done 
in the space 

:= {(u\n'^) G Hi(M" \ F) x \ F) : ([yu"], [yni) G Y, x Y,}, 


which incorporates the only homogeneous essential transmission conditions of (3.4). 










































Proposition 3.2 (Equivalent variational formulation). Problem (3.4) is equivalent to 
finding 


(u^, G H/i s.t. i?((u^, n^), (w, ta); s) = ^((w, w); s) V(w,w)e]HI/, 


(3.7) 


where 


B{{u,v),{w,w)]s) := (cr(u),£(w))K<i\r+ Pss^(u,w)ffid 

+P/(Vn, Vw)^d\r + Pf{s/cf{v, w)^d 
+p/s((7^n, N7”w)r — (7“n,N7’^w)r 

+ (N7+u,7“w)r - (N7"u,7+w)r), 

£((w,m;);s) := pj ((Aq, 7’^w)r - s(NVo, 7”w)r) . 

Proof. The definition of the normal traction for u and the normal derivative for n, plus 
simple algebraic manipulations, show that 


P/(Au, w)^d\r + P/(Vu, Vw)^d\r + (A*u, w)K<i\r + (cr(u), £(w))Kd\r 

= ([t(u)], 7 “w)r + (t+(u), [7w])r + p/( 9 “u, [ 7 '«;])r + P/([M, 7 ^w)r 
= (t+(u) + pjsN*7"u, |7w])r + p{dfiv + SN7+U, [ 7 w])r 
+ ([t(u)] -p/sN*[7u],7-w)r + p/([M - sN[7u],7+u;)r 
+ P/S ((N*7“u, 7'*^w)r - (N‘7 ’^u, 7~w)r + (N7“u, 7’^w)r - (N7+U, 7“tn)r) , 


or equivalently 

B((u, v), (w, w);s) + (A*u - pss^u, w)Kd\r + P/(Au - (s/c)^v, 

= (t+(u) + p/sN* 7 "u, [ 7 w])r + p(dfiv + sNq+u, | 7 w])r 

- P/(sN[ 7 u] - |a^u],7+'«;)r + ([t(u)] - p/sN*[ 7 u], 7 'w)r. 


(3.8) 


From here it is clear that a solution of (3.4) satisfies (3.7). Reciprocally, if we have a 
solution of (3.7), testing with smooth functions with compact support in \ T, we can 
easily recover equations (3.4a) and (3.4b). Therefore, by (3.8) it follows that 

-p/(sN[7u'*] - [d^v^\ + Ao,7^'«^)r + ([t(u'^)l - p/sN^dyu'*] + 0o),7"w)r 

+ (t+(u'") + p/sN*7"u'', [7 w])r + p{dfv^ + sNy+u'*, [7w])r = 0, 


for all (w,tn) G Mh. The transmission conditions (3.4c), (3.4d), and (3.4f) follow from the 
simple observation that the map 

H^/\T)xU^/^{T)xYhXYh 
(w, w) I—^ (y+w, 7 -w, lyw], [yw]) 


is surjective. □ 

The third step in the analysis is the proof of well-posedness of the variational problem 


(3.7). Following [29], we define the energy norm 

III Ill PI := (<T(u),£(n))K<7\r + ||sv^u||2^ + p/||Vu||K.i\r + P/c“^l|s'y| 
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We will also write a := Res > 0 (for all s G C+) and a := min{(T, 1}. To shorten some of 
the forthcoming expressions, we will denote: 


II(0,0)IIW 

||(A, 


= (<^(u),£(u))r,^p + II^uIIr, +pj||Vn||Kd\r + P/c H-^ 

= ll0llW + ll‘^lliAr> 

2 I ll,,l|2 


-1/2,r 


+ 


1/2,r- 


Note that the energy norm and the hrst of the above norms are related by 

|s| 




<i(u,n)|p|<^| 




(3.9) 


Finally, the expression independent of h will be used to mean independent of the choice 
of the spaces Yh and W. 


Proposition 3.3 (Well-posedness). Problem (3.7) is uniquely solvable for any (Ao,0o) ^ 
j^-i/2(r) X i7i/2(r) and s G C+. Moreover, there exists C > 0, independent of h, such 
that 


(u^^;'*) 


||s| < C'^T^IKAq, S0o)||-l/2,l/2,r- 


Proof. A simple computation shows that 


and that 


Re (sB((u,ti),(u,'ti);s)) = ff|||(u, !i)|||f,|, 
u,ti);s)| < C^)/||(Ao,#o)||-i/ 2 ,i/ 2 ,rl|(u,ti)||i,[[j\r, 


(3.10) 

(3.11) 

(3.12) 


which proves well-posedness of (3.7) by the Lax-Milgram lemma. The estimate (3.10) is 
a direct consequence of (3.11) and (3.12), using (3.9) to relate the norms. □ 

The hnal step wraps up the analysis by collecting information from the previous results. 


Corollary 3.4. Equations (3.1) are uniquely solvable for all s G C+ and any choice of the 
closed spaces Yh and Yh- Moreover, if are defined using (3.3) from the solution 

of (3.1), the following bounds hold with C > 0 independent of h: 

< C 




11 ,: 


CT^CT 


II (Ao, s0o)|| -i/2,i/2,r, 


ll(0Ei </’/)lli/ 2 ,r < C'^^IKAo, S0o)||-i/2,i/2,r- 


(3.13a) 

(3.13b) 


Proof. Propositions 3.1 and 3.2 relate the discrete integral system (3.1) to the variational 
problem (3.7), which is shown to be uniquely solvable in Proposition 3.3 The estimate 
(3.13a) follows from (3.10) and (3.9). Finally, the bound (3.13b) follows from (3.13a) and 
(T^ □ 


We end this section by noting that Corollary 3^ implies the unique solvability of 
the semidiscrete equations that are obtained by taking the inverse Laplace transform of 
(3.1). They can also be translated into a time-domain estimate that bounds norms of the 
solution in terms of bounds for the data. 
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3.2 The effect of Galerkin semidiscretization 


In this section we analyze the effect of space semidiscretization, that is, we estimate the 


difference between the solution of (2.4) and (3.1). The analysis follows a very similar 


pattern to the one displayed in Section 


3.1[ We start by writing the error equations: 


■#■/)£ YJxiX 


(3.14) 


We will develop the analysis in terms of the variables := — m and := — v, from 

which the error of the boundary unknowns can be recovered: 


(</>l - ips, ‘t’’} - ih) = (he'*], -|7e'‘l 


(3.15) 


The potential representation for (e™,e^) is obtained by subtracting (2.5) from (3.3) 


e' = -p^sS(s)N‘(.),5-.),,)-D(s)(</.|-</.ri, 

e* = sS(s/c)N(0|-<Ae) + D(s/c)«-■(>/). 


(3.16a) 

(3.16b) 


The proofs of the following results are quite similar to those of Propositions 3. HE ^ and 
|3.3[ We will only point out the main differences. 

Proposition 3.5. The error potentials := — m and e^ := — v satisfy: 


A^e'^-pssV = 0 

in Mf' \ T, 

(3.17a) 

Ae'^-(s/c)V = 0 

in \ T, 

(3.17b) 

sNM-M = 0, 


(3.17c) 

|t(e'*)]= 0, 


(3.17d) 

(he'"], [qei) + (0s, -0/) e Yh 

xYh, 

(3.17e) 

(sNy+e^ + dfe’^, t+(e'^) + pjsN^ye^) G Y° 

xYh- 

(3.17f) 


Reciprocally, given a solution of (3.17), the quantities defined in (3.16) satisfy (|3.14[) 


Proposition 3.6. Problem (3.17) is equivalent to the variational problem: find (e^,e^) G 
H^(M'^ \ T) X \ T) such that 


e"]+0s,[7e1-0/)eY,xW, 
B{{e'^,e^), (w,'u;);s) =0 


V(w, tc) G EI/j, 


(3.18a) 

(3.18b) 


We note that, in comparison with (3.7), problem (3.18) has homogeneous right-hand 
side but incorporates a side restriction (3.18a). This compares with how the conditions 
(3.4c)-(3.4d) have become homogeneous in (3.17c)-(3.17d), while the homogeneous con¬ 
dition (3.4e) is now non-homogeneous (3.17e). 

Proposition 3.7. Problem (3.18) is uniquely solvable for any (0s, 0/) ^ ]HI^/^(r) and 
s G C+. Moreover, there exists C > 0 independent of h such that 


(e^e'^) 


lid < C*—ll( 0 s, 0/)||i/2,r- 


(3.19) 
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Proof. Using the definition of the bilinear form B (see Proposition |3.2[ ) and (3.9), we can 
easily bound 


\B{{u,v),{w,w);s)\ < |||(u,'i;)|||p||||(w,'«;)|||p| +(7151 ||(u,r;)||^K+r||(w,'«;)||i,Kd\r 


(3.20) 


<(7—||(u,t;)||i_K+rlll(w,w)||||,|. 

Take now (w, w) G H^(M'^ \ T) x i7^(M‘^ \ P) such that 

hwl = 0s, tfwj = -(j)f, ||(w,'u;)||i,Kd\r < <7||(0s,0/)||i/2,r- 

By the energy identity (3.11), the fact that (e^ + w, + tc) G H/j, and ( |3.20 ), it follows 
that 


(3.21) 


(e'" + w, e'" + w) 


II PI 


< ^|5((e'^ + w,e'^ + w;),(e'* + w,e'* + w;);s) 


a 


= —\B{{w,w),{e^ + w,e^ + w);s) 


a 


< (7^||(w,w;)||i,; 

aa 


e +w,e +w)|||p|- 


Therefore, using (3.9) 


(e'^,e^)||||,| < (7^||(w,w)||i,K+r, 


and the result follows from (3.21). For readers who are acquainted with this kind of 
Laplace-domain estimates, let us clarify that the use of the optimal |s|-dependent lifting 
of Bamberger-HaDuong [U Lemma 1] (see also |35l Proposition 2.5.1]) instead of the plain 
lifting used in (3.21) does not improve the estimate. This is principally due to the s factor 
in the boundary terms of the bilinear form B. □ 


Corollary 3.8. Let (0e,0/) and (0s)0/) respective solutions of (2.4) and (3.1). 

Let then (u,n) and (u^,n^) be defined through (2.5) and (3.3) respectively. Then there 
exists (7 > 0 independent of h such that 


(u'*-u,'y'^-'i;)||i,Rd\r < (7^||(0s, 0 /)||i/ 2 ,r, 


aa^ 


ll(0s ~ 0S)0/~ 0/)lli/2,r < (7^^II (0s, 0/) II 1 / 2 ,r- 


Proof. The result is a direct consequence of Propositions 3.5, 3.6, and 3.7 


□ 


Using the results obtained in the previous two subsections it is possible to establish 
error estimates in the time-domain. Data will be taken in the Sobolev spaces 

y^fc(^±i/ 2 ^r)) _ ^ C^-i(M;i7±U2^r)) : ^ = 0 in (-cx), 0), G L\R- 

for A; > 1. A straightforward application of the inversion theorem of the Laplace transform 
[m Theorem 7.1] (see also [35l Proposition 3.2.2]) starting with the bounds of Corollary 
3. 4| yields the following: 
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Corollary 3.9. If the data of the problem satisfy Aq G W^{H 0o G 

then (0s, 4>f) and (u^, v^) are continuous causal functions of time and for all t > 0 

||(0s,0/)(t)||i/2,r < ^^max{l,t‘^} ||^3(Ao, 0o)(r)||-i/2,i/2,r 

||(u'*,'i;^)(t)||^K‘^\r < ll^3(Ao,0o)(^)||-i/2,i/2,r c?r, 

where Di and D 2 depend only on T and 


in/m 


In a similar fashion, a combined application of HB Theorem 7.1] and Corollary |3.8 
provides the following estimate for the errors of semidiscretization in time. Note that we 
are allowed to insert the best approximation operators in the right-hand side of the bound 


of Corollary |3.10| because the error produced by trying to compute the exact solution and 
the difference of the exact solution with its best approximation is the same. 


Corollary 3.10. If the exact solution of (|2.4|) satisfies 0s G lT^(H^/^(r)) and (j)f G 
then (e'*,e'*) := (u-u^,n - “g C(M, \ F) x H\R^\T)) and for 

all t > 0 we have the bound 


Dt^ /■* 

||(e'*,e'*)(t)||i_K<^\r < T^maxll,^^} / ||p4(0s - n,,0s, 0/- nft0/)(r)||i/2,r dr, 

^ + J- Jo 

where II/i and fl/j are the best approximation operators inY^ and Y^, and D depends only 
on F. 


3.3 A fully discrete method 


Full discretization with BDF2-CQ. A fully discrete method can be obtained by 
using any of the many Convolution Quadrature schemes. The reader is referred to mu 
122] for the algorithmic description of multistep and multistage CQ schemes. We next give 
an estimate for the BDF2-based CQ method, based on the stability bound in the Laplace- 
domain obtained in Proposition 3.3 and [35l Proposition 4.6.1] (a slight rehnement of one 
of the main convergence theorems in |3nj). 


Proposition 3.11. Let i = 6 and (Aq, 0o) be causal problem data such that Aq G Wf^{H ^'^^(P)) 
and 00 G lF^+^(iFV2(r)), Then 


\\[vL\v'')(t) - < Dk\1 + C) r||(A<'',rf+‘>)(r)||_v 2 ,.Ar*. 

Jo 

It is important to note that the high-order regularity £ = 6 is only required to achieve 
optimal convergence of order k^. For problem data with regularity as low as £ = 3, 
reduced convergence of order is achieved (see [35]). 
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4 General linear elastic materials: BEM-FEM 


Going back to the system of equations (2.3), an alternate approach aiming for a hnite 


element solution of the elastic waveheld and a boundary element solution of the acoustic 
waveheld is to use a direct boundary integral representation of the acoustic wave while 
keeping the partial differential equation for the elastic displacement in variational form. 
This approach is particularly well suited for the case of variable elastic density and Lame 
coefficients, and also for heterogeneous anisotropic materials. In the following we assume 
that the stress is given by a linear law cr = C(x) e, where for each x G G_, C(x) is a 
linear operator that transforms symmetric matrices into symmetric matrices and satishes 
e : C(x)£ > CqS : e for some positive constant Go, for every symmetric matrix e and for 
almost every x G G_. We also assume that the components of the tensor C are L°°{VL_) 
functions and that the solid density ps ^ is strictly positive. The Navier-Lame 

operator is now given by A*u = V ■ (C£(u)) and the traction operator t is redehned 
accordingly as well. 

The derivation employs standard arguments of boundary integral equations and is 
presented with careful detail in [27] , with the resulting equivalent system being 


— A*u = 0 

t~(u) + p/sN*0 = — p/sN*0o 
V(s/c)A+ (|l-K(s/c))0=O 
(—|l + K*(s/c)) A + W(s/c)0 — sNyu = Aq 


in f2_, 
on T, 
on T, 
on T. 


(4,1a) 
(4.1b) 
(4.1c) 
(4. Id) 


For notational convenience, we introduce the interior elastodynamic bilinear form in the 
Laplace-domain 

a(u, w;s) ;= (cr(u), £(w))n_ + s^(p/u, w)q_. 


so that the variational formulation of (4.1a)-(4.1b) reads 


a(u, w; s) + s{pf{(j) + 0o)i 7 w ■ ^')r = 0 Vw G H^(G_). 

We note that the operator Nyw = yw ■ u appears in this weak formulation, while N* will 
not be used any longer in this section. Since the language of this section is less heavy on 
the side of operators, we will keep the explicit form of the combined operator Ny as a 
trace operator dotted with the normal vector held. 


4.1 Galerkin semidiscretization in space 


Just as in Section 3.1, the solvability and stablity of (4.1) are studied simultaneously. In 


order to do so, we dehne the closed subspaces 

V;, C Xf, C Yh C 

The following result establishes the connection between the discrete counterpart of prob¬ 


lem (4.1) and a non-standard transmission problem. Note that the ‘Finite Element’ form 


is a discretization of the interior Navier-Lame equation, and therefore, the elastic op¬ 
erator has been discretized, as opposed to what happens with the ‘Boundary Element’ 
counterpart, where only transmission conditions are discretized. 
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(4.2a) 

(4.2b) 

(4.2c) 

(4.3) 


Proposition 4.1 (Transmission problem for Galerkin equations), //(u^, 0^, A^) G V/^x 
Yh X Xh satisfies the Galerkin equations 

w; s) + s{pf{(j)'^ + (f)o),jw ■ iy)r = 0 Vw G V,* 

-s-fu^ ■ + W{s/c)(j)^ + (-|I + K\s/c)) \^-\o e 

(|l-K(s/c))0" + V(s/c)A" G X°, 

and 

v'^ := D(s/c)0'* - S(s/c)A^ 

then the pair (u^,n^) G V/j x if^(M‘^ \ T) satisfies the transmission problem 
a(u^, w;s) + s(p/(-[7^'^l + 0o),7w ■ j/)r = 0 Vw G V/,, 

(s/c)V =0 
[7n'*l Gn, 

Id^l gx,, 

570^* ■ jy + dfiv’^ + Ao G 

GX, 


m \ r, 


'-h' 


(4.4a) 

(4.4b) 

(4.4c) 

(4.4d) 

(4.4e) 

(4.4f) 


Conversely, given a solution of (4.4), the triplet 

(u^./.A") := e V» X n X ATi 


(4.5) 


satisfies (4.2). 


Proof. Equations (4.4b), (4.4c), and (4.4d) are simple consequences of the definition of 
and the jump relations of the double and single layer potentials. Moreover, using the 
dehnition of and the well known identities 

a;S(s) = + Kfis ), 7 'D(s) = -^I + K(s), 


it is easy to verify that (4.4e) and (4.4f) are just restatements of (4.2b) and (4.2c). 

To prove the converse, note that (4.4b) and the dehnition of in (4.5) imply the 

integral representation (4.3). Then (4.2b) is equivalent to (4.4e) and (4.2c) is equivalent 
to ( |Of| ). □ 

Proposition 4.2 (Equivalent variational formulation). Let 

V, := {n G H\M.^ \ T) : [yn] G Y, , yn G X°}. 


The problem (4.4) is equivalent to finding (u^,n") G V/j x 14 such that 


where 


and 


A{{u’^,v’^),{w,w)-,s) = f{{w,w)]s) V(w,u;) G x 14, 

A ((u, n), (w, w); s) := (<t(u), £(w))^_ + (p^u, w)^_ 

+ Pf (Vn, Vu;)Kd\r + pAs/cf {y, u;)Kd\r 

+ P/s(7u ■ 1^, - p/s(|7n], 7w ■ jy)r, 

/ ((w, w); s) := -pfs{(j)o, yw ■ jy)r - p/(Ao, [7w])r- 


(4.6) 
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Proof. Let be a solution pair for (4.4). Then, for all w G 14, 


w)r - {dfv^,'y^w)r - w)r 

= (Vn^, Vw)]Rd\r + {s/cf{v^, w)^d, 

after applying (4.4b) and ( 4.4f| ). Therefore, testing (4.4e) with [yw] for w E Vh, and 
substituting the above, it follows that 

{s/cY{v^,w)^d + (yv'^,Vw)^d\r + [7H)r = -(Ao, [7H)r Vw G I 4 . (4.7) 


However, the pair of equations (4.4a) and (4.7) are equivalent to (4.6). 

To prove the converse statement, note that we need to show that a solution of (4.7) sat- 
ishes (4.4b), (4.4d), and (4.4e). Equation (4.7) applied to a general compactly supported 
w G C“(M“\r) is the distributional form of (4.4b). Therefore, (4.7) (after integration by 
parts) implies 

{Of o’", 7"w)r - 7+w)r + (s7U^ ■ + Xo, [7^)r = 0 Vw G 14, 

which, after some simple algebra, is shown to be equivalent to 

{dfv'^ + S7u^ ■ + Ao, [7^)r + {ld^v^j,'y~w)r = 0 Vw G 14- (4.8) 

However, the operator Vh 3 w \— > (| 7 tn], 7 “tn) E Yh x is surjective, and therefore 
(|4.8|) is equivalent to (|4.4d[) and (|4.4e|). □ 


For the analysis of ( 4.6[ ), we need to redehne the energy 

lll(u,'y)|||p| := (<T(u),£(u))n_ + ||sy^u||^_ + P/llVu||K.\r +P/ll(s/c)u| 


norm 
|2 


due to the fact that the elastic held is not handled with a potential representation and, 
therefore, it does not extend to the other side of the interface. Note that ||| ■ |||i is equivalent 
to the H^(ff_) X i7^(M'^ \ T) norm and that, similarly to (3.9), 


ctIIKu.iOIIIi < III(«.■<')Illi.i < yIII(u.I')111- 


(4,9) 


Proposition 4.3 (Well-posedness). Problem (4.6) is uniquely solvable for any (0o, Aq) G 
7/^/^(r) X and s E C+. Moreover, there exist Ci, C 2 > 0, independent of h, such 

that 


|||(u^,u^)|||i + || 0 ^||i/ 2 ,r < 41 —^||(s0o, Ao)||i/ 2 -i/ 2 ,r, (4-10) 

era"' 

||A II - 1 / 2 ,r < 42 ^^ 3/2 II Ao) II 1 / 2 - 1 / 2 ,r- (4-11) 

Proof. It is straightforward to verify that 

Re(s7l((u,u),(u,u);s)) = a|||(u, u)|||f,|, (4.12) 
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and 


\f{{w,w);s) I < CIKs^o, Ao)||i/ 2 ,-i/ 2 ,r|||(w,'w;) 


where the constant depends only on pf and F. Hence, by (4.9) and the Lax-Milgram 
lemma, we have nniqne solvability of (4.6) and the following bonnd in the energy norm: 

(4.13) 


(u ,n )||||s| < C*—||(s 0 o, Ao)||i/2,-i/2,r- 


The estimate (4.10) can be easily derived from (4.13) and (4.9) and the fact that 0" = 
— Finally, recalling that and nsing [2^ Lemma 15], namely if Av — 

s^v = 0 in an open set O with Lipschitz bonndary, then 


/I 

\\duv\\-i/2,do < C ( ^ j (lls^'llo + llVnllo), 


it can be shown that (4.11) follows from (4.10). 


(4.14) 

□ 


4.2 Semidiscretization error 


We now stndy the difference between the solntions to the exact problem and their hnite 
dimensional approximations. It is important to stress that — u ^ V/i, and therefore 
we will not be able to write a transmission problem for the error — u in the style of 
(4.4). Instead, we will work with the difference 


e := u 


PftU, 


where P/i : H^(r2_) —>■ V/^ is an elliptic projection that will be dehned below. We hrst 
need to introdnce the hnite dimensional space of rigid motions 

M := {m G H^(H_) : (cr(m), £(m))Q_ = O} . 

From now on we will assnme that M C V/^. The operator P/^ is given by the solntion of 
the problem 


(cr(P;ju),£(w))n_ = (cr(u),£(w))n_ Vw G V,,, (4.15a) 

(PfeU, m)n_ = (u, m)Q_ Vm G M. (4.15b) 

Using Korn’s ineqnality it is easy to show that P/^ is well dehned and that the approxi¬ 
mation error ||u — P/iu||iq_ is eqnivalent to the H^(r2_)-best approximation on V/j. In 
order to shorten notation, we will write := P/^u — u. 

The triplet (e^, 0^, A^) G V/j x Yh X Xh satishes the following error eqnations: 

a(e'', w;s)n_ + (psr'*, w)^_p/s((0'* - 0),7W • i/)r = 0 Vw G V/^ (4.16a) 

+ y'^)-u + W(s/c)(0'^ - 0) - (|I - K*(s/c)) (A'* - A) G (4.16b) 

ill - K{s/c)) (0" - 0) + V(s/c)(A" - A) G x;;. (4.16c) 

For this system there is a corresponding non standard transmission problem 
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Proposition 4.4. //(e",A",0") satisfies (4.16) and we define 

e^ := D(s/c)(0'* - 0) - S(s/c)(A'^ - A), 
then the pair then (e^, e^) is a solution of the transmission problem 


+ {s/cfe^ 

-hll A 


l'ye‘ 




)Je'- 


= - (psr^, w)^_ VwGV/i, 

(4.17a) 

= 0 m \ T, 

(4.17b) 

G^, 

(4.17c) 

eXh, 

(4.17d) 

elk, 

(4.17e) 

eX4 

(4.17f) 


Conversely, i/(e^,e^) is a solution of (4.17) then 


(e^0^A") ;=(e^0-[7e"],A-[a.e"l), 


solve (4.16). 


Proof. Starting with a solution of (4.16), we see that equation (4.17b) is a consequence of 
the dehnition of e^, while (4.17a) follows readily from (4.16a), noting that ([ye^], = 

(0 — 0^, A — A^). The equations (4.17c) and (4.17d) can also be verihed from the last 
observation, since Yh x Xh 3 XV = (0 — [qe^], A — [9i,e^]). Finally, using 

a;(S(s)A) = (il + K‘(s))A, 7-(D(s)4,) = (-11 + K(s)) 4 ,, 

we see that (4.16b) and (4.16c) imply (4.17e) and ( |4.17f| ). 

The proof of the converse statement is very similar. □ 

Proposition 4.5. The system (|4.17) is equivalent to the variational problem of finding 
(e'^, eV e X \ T) such that 

(^-e^ [qel - 0) e X n, 

Al((e'",e'^), (w,r(;);s) = 6((w,r(;);s) V(w,w) G V/, x 14, 
where the bilinear form A is defined in the statement of Proposition o,nd 
b ((w, w); s) := pf{X, 7“w)r + sp/( 7 r'* ■ ly, [7w])r - (psr'*, w)^_ 


(4.18a) 

(4.18b) 


Proof. The proof is very similar to the one of Proposition 4.2 Details are omitted. □ 


Proposition 4.6. Problem (4.18) is uniquely solvable for any (n, 0, A) G H^(r2_) x 
i7^/^(r) X i7“^/^(r) and s G C+. Moreover, there exist constants C\,C 2 > 0 independent 
of h such that 


\ e^)|||i + \\fi - 0^111/2,r (^Il(s<(>, A) II1/2,-1/2,r + 


Isr'^lko- + ||sk'^|k_ 


(4.19) 


\\X — A^||i/2,r <C 2 ^ j.3/2 (iK-^ II i/2,-i/2,r + ll-sr^lkn. + ||s^r^||n_j. 


(4.20) 
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Proof. The existence and uniqueness of the solution to (4.18) is proven in a way analogous 


to that used in Proposition 4.3 We will next prove a bound in the energy norm 


'^,e^)||||s| < Cl —^11 (s0, A)II 1 / 2 - 1 / 2 ,r + ||■sr^||l,Q_ + ||s 


2„h|| 

r hi- 


(4.21) 


The estimate (4.19) follows from (4.21) and (4.9). In order to get to (4.20) we make use 
of (4.21), the fact that A — = [c?ye^], and (4.14). 


To prove (4.21) we proceed as follows. We first obtain an upper bound for the bilinear 
form 


IAl((u, n), (w, w); s)I < C^III(u, v )||| 11 ||(w, w )||||,|, 


by the same argument that was used in Proposition 3.7 Also 


C 


\b{{w,w)]s) I < - (||A||_i/ 2 ,r + llsr^lko- + ||s^r'^lb_) |||(w,m;)||||,|. 


(4.22) 


(4.23) 


The constants in (4.22) and ( |4.23 ) depend only on the geometry. Now, for 0 G 
pick a lifting w^f, G \ P) such that = 0, 'y~w^ = 0, and 


ll'U^0||i,Rd\r < C||())||i/2, 


r- 


(4.24) 


Since (e^,e^ + w/,) G V/^ x 14 we can use (4.12), (4.18b), (4.22), and (4.23) (i.e., coercivity, 
the variational equation, and boundedness of the bilinear form and right-hand side) to 
estimate 


(e'^,e^ + W0)p,l <^|Al((e'*,e^-F'«;,/),(e'*,e'^ + w;^);s) | 

= — \b ((e'*,e^ + tc/,);s) -h ^ ((0, w,^), (e^, e'* -1-w^); s) 


a 


<C^|||(e4e" + ^^)||||,| 


aa 


(klll'*^0lli,iR‘^\r + ||'^||-i/2,r + ||■sr^||l,Q_ -|- ||s^r^|4_) 


This bound, together with 


1(0,0;^) III PI < —||s0||i/2,r 


see (4.9) and (4.24)) prove (4.21) 


□ 


4.3 Estimates in the time-domain 

Using the bounds for the error operators derived in the previous section, we can prove 
explicit time-domain estimates. Just like in the BEM/BEM case, we can use [HI Theorem 
7.1] and combine it with the Laplace-domain estimates from Propositions 4.3 and 4.6 to 
obtain the following results. 
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Corollary 4. 7. Consider causal problem data <1)0 € and Xq G Wl{H ^/^(F)). 

Then u^, are continuous causal functions of time and for all t > 0 : 




+ || 0 ^Wlli/ 2 ,r <Di max{l,t^} 


t + 1 


||A'*(t)||_i/ 2 ,r <£> 2 max{l,t^/^} 


ty/i 
y/t + 1 


||’P3(0O5-^o)('r)||i/ 2 ,-i/ 2 ,r dr, 
rt 

||’P3(0O5 Ao)(r)||i/ 2 ,_i/ 2 ,r dr. 


where Di and D 2 depend only on F. 


To abbreviate the following statements, we will use the following shorthand for ap¬ 
proximation errors 


«/.(«) - j‘ {md - ni'^'-)(r)||i/ 2 ,r + ||P3(A - nf A'*)(r)|U,/ 2 ,r) dr 
+ [ (ll^3(u-P/iu'^)(r)||i,n_ + ||p3(u-P/iu'*)(r)||n_) dr, 


where Yh and : H' 

the elliptic elastic projection onto Yh dehned in (4.15). 


^/2(r) ^ Xh are orthogonal projections and Ph is 


Corollary 4.8. If the solution triplet satisfies 


(u,0,A) G hF^(H^(fi_)) X W^{H^/\r)) X 

then (e^,e^) G C(M, H^(r2_) x \ F)) is causal and we have constants Di and D 2 

depending only on F such that for t > 0 

|||(e^,e'*)(f)|||i + ||( 0 - 0 '^)(f)||i/ 2 ,r < T>i max{l, ah{t), 


+3/2 

||(A - A'^)(t)||_i/2,r < T>2max{l,t^/^}—== Ohit). 

y/t + 1 

Full discretization with BDF2-CQ The purely boundary integral formulation treated 
in the first part of this paper lent itself naturally to a full discretization using one of the 
many Convolution Quadrature schemes for the time evolution. For the current varia¬ 
tional/boundary integral formulation it would seem that an independent treatment with 
traditional time-stepping for the Finite Element part and Convolution Quadrature for the 
discretized boundary integral equations would be the best way to proceed, and for our 
computational implementation we will proceed in this fashion. 

However, it turns out that the separate application of time stepping and CQ to different 
parts of the system is equivalent to the application of CQ globally, as long as the time 
stepping method used for the FEM part coincides with the one giving rise to the CQ 
family used for the implementation (see [29l Proposition 12], [23]). This observation will 
allow us to analyze the fully discrete method as if the whole discretization were done with 
CQ. 
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We present results for the coupled schemes based on BDF2. In the following section 
we will show numerical experiments for BDF2-CQ and Trapezoidal Rule-CQ. (We note 
that the analysis of Trapezoidal Rule CQ was done by Lehel Banjai in [2], although it does 
not give explicit behaviour of bounds with respect to t.) We will use to denote 

the fully discrete approximation of (u, v) using a CQ method with constant time-step k. 
In parallel to the corresponding result in Section]^ (Proposition |3.lf| , the next estimate 
follows from the Laplace-domain estimates in Proposition |4.6| and an application of [3^ 
Proposition 4.6.1]. 

Proposition 4.9. Let i = 6 and (0o, -^o) be causal problem data such that 

(00, Ao) e X Wi{H-^/\T)). 


Then, for t > 0, it holds that 

|||(u'*,n^)(f) - (u^,'y;^)(t)|||i < T>K^(l + fQ f ||(0o^^^\A[,^^)(r)||i/2,-i/2,rdr. 


It is important to note that the high-order regularity £ = 6 is only required to achieve 
optimal convergence of order k"^. For problem data with regularity as low as i = 3, 
reduced convergence of order is achieved (see [35], Chapter 4]). 


5 Numerical Experiments 

In this section, we show some experiments for fully discrete methods applied to the BEM 
and BEM/FEM formulations we have analyzed. For general ideas of what CQ time- 
discretization means and how it is used, we refer to [3,1221 IS]- Algorithms for BEM/FEM 
applied to acoustic transmission problems are explained in [23] - 


5.1 Boundary integral method 

In order to test the convergence properties of the implementation the following synthetic 
problem was solved in The interior elastic domain will be the unit disk = {x : 
xl + < 1}, and its exterior will be the acoustic domain. If we let 'R(f) be a smooth 

approximation to the Heaviside function, then the elastic causal pressure wave 


u(x,f) ='H(cLf-x-d)sin(3(cLf-x-d))d, d = , cl = 

and the cylindrical acoustic wave 

^(x, t) = C~^ |«ilQ^^(*s|x|) C{'H{t) sin(2f)}| 


_ / 2/iH~A 


solve equations (2.2a) and (2.2b) respectively. Here C is the Laplace transform. In order 


for them to satisfy the entire IBVP (2.2), equations (2.2c) and (2.2d) were used to dehne 


the boundary data ao := and /3o := n* 
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The boundary data was sampled accordingly and the Laplace transformed equivalent 
system (2.4) was discretized in space with deltaBEM (the reader is referred to [101 E] for 
further details on the computational aspects of deltaBEM), which can be considered as a 
Galerkin Vi method with reduced quadrature, while Convolution Quadrature was used 
for time stepping on increasingly hner space/time discretizations with N space points and 
M time steps. The approximated solutions were then sampled in 20 random points on the 
circle of radius r = .7 for the elastic wave and r = 2 for the acoustic wave and compared 
against the exact solutions. The maximum difference in the hnal time 


max2£i 

\v{Xi,tf) - 

max‘f 2 ,i 
. max2£ju(xi,t 



|u(Xi,t/)| 


is used as the error measure. Trapezoidal Rule CQ and BDF2-CQ were both implemented 
and compared. Tables [T] and [^summarize the results, while convergence plots can be seen 
in Figure [2 In the simulations, the values A = 9, /i = 15, ps = 1-5, p/ = 1 and c = \/5 
were used, the hnal time was T = 5. 


N/M 

^h,k 

e.c.r. 

T?V 

e.c.r. 

45/90 

0.8745 

— 

1.1603 

— 

60/120 

0.7131 

1.2265 

0.9862 

1.1766 

90/180 

0.3692 

1.9312 

0.8900 

1.1080 

120/240 

0.2022 

1.8265 

0.4778 

1.8627 

180/360 

0.0806 

2.5079 

0.2407 

1.9854 

240/480 

0.0466 

1.7285 

0.1482 

1.6240 

360/720 

0.0302 

1.5456 

0.0513 

2.8869 


Table 1: Relative errors and estimated convergence rates in the time-domain for the 
BDF2 Convolution Quadrature with lowest order Calerkin discretization (with reduced 
quadrature): N represents the number of space discretization points (elements), M is the 
number of timesteps. The errors are measured at the hnal time T = 5. 


5.2 Coupled boundary-field method 

The previous coupling scheme was implemented using 'P 3 /P 2 Boundary Elements for 
the acoustic wave held and V 3 Finite Elements for the interior elastic waveheld. The 
convergence studies were carried out using the rectangle [1,3] x [1,2] as the elastic domain, 
which was triangulated using Matlab-produced unstructured meshes. Known solutions 
were imposed for the interior and exterior problems; a plane pressure wave on the interior 

u = V’(cLt - X ■ d) d, ^p{t) := 

and a cylindrical acoustic wave on the exterior 

n ={|i/i^^(3|x-Xo|)£{<p(t)}| , <p(t) := H(f)sin(3t) 
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N/M 

^h,k 

e.c.r. 

jpv 

e.c.r. 

45/90 

9.9416 

— 

1.0642 

— 

60/120 

73.3473 

0.135 

0.3961 

2.6864 

90/180 

0.1402 

523.3082 

0.2089 

1.8962 

120/240 

0.0675 

2.0755 

0.1261 

1.6571 

180/360 

0.0484 

1.3955 

0.0522 

2.4138 

240/480 

0.0252 

1.9219 

0.0308 

1.6968 

360/720 

0.0181 

1.3923 

0.0126 

2.4404 

480/960 

0.0099 

1.8212 

0.0066 

1.9054 


Table 2: Relative errors and estimated convergence rates in the time-domain for the 
Trapezoidal Rnle Convolntion Qnadratnre with the same space discretization as in Table 
[TJ N represents the nnmber of space discretization points, M is the nnmber of timesteps. 
The errors are measnred at the hnal time T = 5. 



Fignre 2: Relative errors for the BDF2 and TR implementations of CQ. The maximnm 
difference between the approximate and exact solntions on the sampled points. 


where xq = (1.5,1.5) is the location of the sonrce of the cilyndrical wave, 'H{t) is a smooth 
approximation to the Heaviside fnnction, and A = 2, p = 3 and p = 5. 

These two fnnctions satisfy eqnations (2.3a) and (2.3b). In order to force them to solve 
the problem in qnestion, the bonndary data was mannfactnred using (2.3c) and (2.3b) 
as the dehnitions for (Ao,0o)- The relevant information was sampled from the known 
solution, combined according to (2.3c) and (2.3b) and the resulting pair (Aq, 0o) was then 
fed to the discrete system as boundary data. 

The experiment was run with hxed FEM and BEM grids with h = 0.025 (maximum 
element area 3.5 x 10“^ for the FEM mesh) for a hnal time T = 1.5. Trapezoidal Rule 
time stepping and Trapezoidal Rule Convolution Quadrature was used respectively for 
the Finite Element and Boundary Element domains with doubling number of time steps 
starting at 5 and all the way up to 160. The errors were measured for the hnal time, for 
the hnite element solution in the L^(f2_) norm and in the Fr^(f2_) norm. 

For the acoustic waveheld the discrete solution was postprocessed, sampled and compared 
to the exact solution in 10 random points in the acoustic domain, with the -normalized- 
maximum discrepancy ^ being considered as the error. 

The last two tables and hgures show the convergence studies for the Trapezoidal Rule 
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Figure 3: Convergence studies for the coupled FEM/BEM scheme. elements were used 
for the finite element solution and 'P^IV 2 elements for the boundary element solution. 


M 

JPV 

e.c.r. 


e.c.r. 


e.c.r. 

5 

4.051 E-3 

— 

4.531 E-2 

— 

0.179 E-1 

— 

10 

3.393 E-3 

0.255 

1.204 E-2 

1.912 

7.188 E-2 

1.315 

20 

1.925 E-3 

0.818 

3.364 E-3 

1.839 

2.524 E-2 

1.510 

40 

5.108 E-4 

1.912 

8.800 E-4 

1.934 

7.184 E-3 

1.813 

80 

1.281 E-4 

1.996 

2.223 E-4 

1.981 

1.862 E-3 

1.947 

160 

3.201 E-5 

2.000 

5.592 E-5 

1.995 

4.700 E-4 

1.987 


Table 3: Relative errors and estimated convergence rates in the time-domain for the 
Trapezoidal Rule Convolution Quadrature. M is the number of timesteps. Final time 


T = 1.5. 


BEM/FEM scheme with simultaneous space/time refinement. Polynomial degrees k = 1 
and k = 2 were used. 


6 Conclusions 

We have presented stable fully discrete formulations for wave-structure scattering in the 
time-domain. The formulations are well suited for pure Boundary Element implementa¬ 
tion or Boundary/Finite element coupling in the case of a general elastic scatterer. For 
time discretization a BDF2 Convolution Quadrature scheme was used as a basis for the 
analysis, but, as the numerical experiments show, the Trapezoidal Rule schemes present 
good global convergence properties. 

The analysis presented in this work generalizes easily to the case of several scatterers 
or those with non simply connected geometries. The study of scattering by obstacles with 
more complex physical properties such as piezoelectrics is the object of current research 
and could be of interest in applications such as active noise reduction. 
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Figure 4: Simultaneous space/time refinements. Left: Vi elements were used for the 
finite element solution and Vi/Vq elements for the boundary element solution. Right: V 2 
elements were used for the hnite element solution and V 2 IV 1 elements for the boundary 
element solution. 


M/N 

JPV 

e.c.r. 

-^h,fc,L2 

e.c.r. 

rpu 

e.c.r. 

10/1 

3.412 E-3 

— 

1.382 E-2 

— 

1.555 E-l 

— 

20/2 

1.929 E-3 

0.823 

4.112 E-3 

1.749 

7.467 E-2 

1.058 

40/3 

5.113 E-4 

1.915 

1.088 E-3 

1.917 

3.604 E-2 

1.051 

80/4 

1.281 E-4 

1.996 

2.763 E-4 

1.978 

1.776 E-2 

1.021 

160/5 

3.202 E-5 

2.000 

6.935 E-5 

1.994 

8.843 E-3 

1.006 


Table 4: Relative errors and estimated convergence rates in the time-domain for the 
Trapezoidal Rule Convolution Quadrature with Vi/Vo boundary elements and Vi hnite 
elements, h = 0.52 x 2~^ is the maximum lenght of the triangulation and M is the 
number of timesteps. Final time T = 1.5. 


A purely time-domain analysis, using the theory of evolution equations HD ESI EH 
might shed sharper convergence estimates, with less regularity required for the solution 
and better control of bounds with respect to the time variable. However, the wave- 
structure interaction problem cannot be written as a second order evolution equation for 
an unbounded operator, because of the presence of time derivatives in the transmission 
conditions and other possible avenues have to be explored. 
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